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ABSTRACT 

We present the first implementation of Active Galactic Nuclei (AGN) feedback in the 
form of momentum driven jets in an Adaptive Mesh Refinement (AMR) cosmological 
resimulation of a galaxy cluster. The jets are powered by gas accretion onto Super 
Massive Black Holes (SMBHs) which also grow by mergers. Throughout its formation, 
the cluster experiences different dynamical states: both a morphologically perturbed 
epoch at early times and a relaxed state at late times allowing us to study the different 
modes of BH growth and associated AGN jet feedback. BHs accrete gas efficiently at 
high redshift (z > 2), significantly pre- heating proto-cluster halos. Gas-rich mergers at 
high redshift also fuel strong, episodic jet activity, which transports gas from the proto- 
cluster core to its outer regions. At later times, while the cluster relaxes, the supply of 
cold gas onto the BHs is reduced leading to lower jet activity. Although the cluster is 
still heated by this activity as sound waves propagate from the core to the virial radius, 
the jets inefficiently redistribute gas outwards and a small cooling flow develops, along 
with low-pressure cavities similar to those detected in X-ray observations. Overall, 
our jet implementation of AGN feedback quenches star formation quite efficiently, 
reducing the stellar content of the central cluster galaxy by a factor 3 compared to the 
no AGN case. It also dramatically alters the shape of the gas density profile, bringing 
it in close agreement with the {3 model favoured by observations, producing quite an 
isothermal galaxy cluster for gigayears in the process. However, it still falls short in 
matching the lower than Universal baryon fractions which seem to be commonplace 
in observed galaxy clusters. 

Key words: galaxies: clusters: general - galaxies: active - galaxies: jets - methods: 
numerical 



1 INTRODUCTION 

The well-known over-cooling problem in galaxy formation is 
encountered in numerical simulations for a range of galaxy 
sizes, spanning dwarves to red ellipticals. Stellar feedback 
mechanisms, such as winds from young stars and super- 
novae explosions are potentially good candidates for ex- 
pelling large amounts of cold gas from galaxies. Unfortu- 
nately, these mechanisms are inadequate for the most mas- 
sive galaxies, where the energy liberated by star formation 
activity is too low to unbind material from their gravita- 
tional potential wells. 

Different theories find ways to partially or totally pre- 
vent the cooling catastrophe in the most massive structures: 
galaxy groups and clusters. Propagating heat with a Spitzer 
conductivity from the outskirts to the central parts of a 
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cluster, thermal conduction co uld solve the cooling p roblem 
for the more massive clusters (|Voigt k Fabianll200^ ). How- 
ever, due to the propagation of electrons along magnetic field 
lines, the conduction is essentially anisotropic leading to in- 
sta bilities (such as the Heat flux Buoyancy In stability, HBI, 
see lQuataer'tll2008l : [Parrish k Quataertll2008l ) in the cluster 
core. These can reorient the magnetic field lines and sto p the 
net inflow of heat tow ards the centre (P arrish et al. 20091 : 
iBogdanovic et al]|2009l V Pre- heating of the gas destined to 
fall into a cluster potential well during proto-galactic stages 
has also been proposed as a means to empty ^as on galactic 
scales in these massive halos (|Babul et al.ll2002l ). 

However, in most models of massive galaxy forma- 
tion, AGN play a crucial role in regulating their gas con- 
tent. Red and dead galaxies commonly exhibit the sig- 
natures of SMBHs which are thought to power high ve- 
locity jets into the hot surroundings of the galaxies. Ob- 
servational evidence for strong AGN activity in groups 
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and clusters is plenty ( 


Arnaud et all 19841: Carilli et alJ 


19941: McNamara et alJ 


20011. l2005l: iFabian et all l2002l: 


Birzan et alJ 12004: iForman et alJ l2007l). When sDatiallv 



resolved, this a ctivity often takes the form of radio 
lobes or cavi t ies (iBoehringer et al" 1 19931: 1 Owen et aL I2OOOI: 

■ l2C ■ ■ 



Birzan et aD l20o3: [McNamara et al.1 l2005l: iFabian et al 
iTavlor et"al1 I2OO6I : IPong et al.1 I2OI0I: JPunn e~ 



2OIOI ). or thin extended iets (iBridle et al.l[l994 i). The radio 



lobes or cavities are associated with low AGN activity, i.e. 
the radio mode, and jets with the quasar modes. 

AGN feedback is invoked to efficiently suppress star 
formation in the most massive galaxies either by ejecting 
gas from their Interstellar Medium (ISM) into the Intra- 
Cluster Medium (ICM), or by pr eventing the IGM gas 
from collapsing into galactic discs (|Binnev Sz Taborl Il995l : 
iRephaeh k Silklll995l ). The powerful ejection of gas by AGN 
is also supposed to suppress the formation of cool cores in 
some fraction o f galaxy clusters, a nd in turn to make the 
ICM turbulent ([Dubois et al.ll2009l ). 

Many of the previous numerical studies that attempt 
to simulate the growth of BHs and their associated AGN 
feedback have been performed using t he Smoothed Par- 
ticle Hydrodynamics (SPH) technique (ISiiacki et al.l l2007l : 



Pi Matteo et al .11 20081 : ISooth Schavell2009l ). As shown by 
Mitchell et al.l (|2009h . SPH codes suffer from underestimat- 
ing the true entropjo profile in cluster cores because of their 
trouble resolving Kelvin -Helmoltz instabilit ies in regions of 
strong density contrast (jAgertz et al. I I2OO7I ). As it is gener- 
ally assumed that BHs accrete gas at a Bondi rate which 
is related to the local entropy Mbh oc K~^^^, a poor es- 
timate of this entropy leads to an incorrect calculation of 
the accretion rate onto a BH and, thus, its energy release. 
To circumvent this issue, we use a sink particle approach 
to follow the growth and AGN feedback of BHs within an 
Adaptive Mesh Refinement (AMR) code. 

Moreover, in previous cosmological simulations, AGN 
feedback is modeled with a thermal input of energy. Mean- 
while, theoretical work and observations suggest that AGN 
feedback is mostly mechanical, not thermal. Numerous 
numerical simulations have implemented and tested the for- 
mation and propagation of AGN jets on cluster scales and 
their impact on t he IC M using either idea li zed simulations 
(IChurazov et al.1 I2OOII: iQuilis et al.l I2OOII: iRevnolds et al.1 



200 ll: i Basson & Alexander 
Ruszkowski et al. 2004b; 



'2003; 



Cattaneo 



Gaibler et al 



eyssie: 



a 



Omma e t al 
Ve rnaleo Revnold; 
[2007: ISimionescu et al 



2004 



2006 



2009 



20091 : lO'Nein JonesI I2OIC), or cosmo- 
logica l simulations (iHeinz et al.1 I2OO6I : iMorsonv et al.l 
I2OIOI I. but none so far have followed the BH growth 
self-consistently (i.e. resolving both the BH growth and 
jet- AGN feedback in cosmological simulation over a Hubble 
time). As the AGN feedback is tightly linked to its BH 
growth history, it is of crucial importance to model both 
the BH evolution through time and its jet-energy release. 
In this paper, we propose to bridge this gap by performing 
the first cosmological simulation including a self-consistent 
treatment of BH evolution and its associated AGN jet 
energy release. 



^ Here and later, the entropy is defined as K = Tp where 
T is the gas temperature and p is the gas mass density 



The paper is organized as follows. In section [21 we de- 
scribe the physical ingredients of our simulation. We start by 
presenting our cooling and star formation prescriptions and 
then we introduce the scheme for following the formation 
and mergers of BHs, gas accretion onto them, and their en- 
ergy release in the form of jets. In section (3] we describe our 
initial condition set-up and our simulation run. In section [4l 
we show how the BH growth is linked to the galaxy clus- 
ter formation history, and what drives the different modes 
of AGN feedback (radio versus quasar). In section O we 
demonstrate that this type of anisotropic mechanical AGN 
feedback is able to suppress the cooling catastrophe occur- 
ing in the galaxy cluster. Finally, in section (6] we discuss 
our results. 



2 MODELING THE PHYSICS OF GALAXY 
FORMATION 

2.1 Modeling star formation 

Gas in our simulation is allowed to ra diate energy by atomic 
collis ions in a H/He primordial gas ([Sutherland Dopital 
Il993h so that it ca n collapse i nto dark matter potential wells 
to form galaxies (jSilkl [l977h . To model reionization from 
z — 8.5, heating fr om a UV background is f ollowed with the 
prescriptions from iHaardt k, Madaul (|l996h . Star formation 
occurs in high density regions p > po (po — O.lH.cm"^). 
When the density threshold is surpassed, a random Poisson 
process spawns star cluster particles according to a Schmidt 
law = ep/tff, where tff is the gas free-fall time and e is the 
star formation efficiency, taken to be e = 0.02 i n order to 
repro duce the observational s urface density laws (|Kennicu"ttl 
'1998') . The reader can consult "Ra sera k, TevssieJ (|2006l l and 
Duboi s k T e vssier. (2008b) for more information on the star 
formation implementation. 

In this work, we do not model supernova feedback. Sev- 
eral authors have argued that super no vae can only have a dy- 
namical impact on low mass galaxies (jSpringel k HernquistI 
.2003a : iPubois k Tevssie r 2008b). Thus, as a first order ap- 
proximation, we assume that supernova feedback has very 
little effect on the growth of BHs in massive galaxies as 
they alone appear incapable of removing a substantial frac- 
tion of gas from the ISM. Whilst this simplification allows 
us to properly isolate the effect of AGN feedback on the sur- 
rounding gas from any other galactic feedback mechanism, 
it does not allow us to study the role of metal enrichment on 
cooling in the centre of massive halos. However, we will see 
that the main features of a cooling catastrophe are already 
captured with zero metallicity cooling. 

Finally it is possible to view the modification of the 
temperature at high density p > po by a polytropic equa- 
tion of state (EoS) that we introduce for numerical reasons 
in section 12.31 as a way to take into account the thermal 
effect of the heating of the ISM by supernovae. As a mat- 
ter o f fact, a similar EoS approa ch is used by other authors 
(e.g. ISpringel k HernquistI l2003l ) as a simple model for the 
unresolved multiphase structure of the ISM in cosmologi- 
cal simulations. More specifically, the minimum tempera- 
ture in dense regions becomes Tmin = To(p/ po)^~^ , with 
To = lO'^K, and n — 4/3 which leads to a constant Jeans 
mass Mj = 1.3 10^ M©. Such a value of the polytropic index 
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n roughly compares with the complex functional form of the 
EoS obtained by analyti cal considerations on the mu ltiphase 
structure of the ISM in ISpringel HernauistI (|2003l ). 



2.2 SMBHs as sink particles 

Sink particles were first introduced bv lBate et alJ (|l995l ) in a 
SPH code. Sinks are massive particles that capture gas par- 
ticles in their surroundings. They mimic the formation of un- 
resolved compact objects, e.g. proto-stellar cores in the ISM, 
black holes in the ISM, central super-massive black holes in 
galaxies, etc. Due to the very Lagrangian nature of the sink 
particle technique, it h ad been extensively an d exclusively 
used in SPH codes until lKrumholz et al.1 (l2004l) extended its 
use to grid codes. The ve rsion in RAMSES (iTevssierlliooi ) IS 
strongly inspired by the iKrumholz et al.l ( 20041 ) numerical 
implementation. 

Sink particles are created in regions where the Jeans cri- 
terion is violated, i.e. in regions where the maximum level 
of refinement is reached and where the gas density is large 
enough to potentially produce a numerical instability, in 
other words where: 



Ax 



> Aj 




(1) 



Here Ax is the size of the smallest cell, Aj the Jeans length, 
Cs the sound s peed and p the gas density. According to 
iTruelove et aD |l997h . the numerical stability of a gravita- 
tionally bound object is ensured if it is resolved with at least 
4 cells. With a mixed composition of matter (dark matter, 
gas, stars), Jeans stability is not trivial anymore, but we 
can reasonably assume that gas is the dominant source of 
gravitational potential inside dense collapsed objects, like 
galaxies, in our case. 

For numerical stability, each time that the Jeans crite- 
rion is violated we should spawn a sink particle with a mass 
corresponding to the depleted mass. However, in cosmolog- 
ical simulations this leads to excessively large sink masses. 
The reason is that the gas is concentrated in structures 
(galaxies) that are poorly resolved (kpc scale). As a result 
an entire galactic disk can be defined by only a few Jeans- 
violating cells leading to massive sinks. To form sufficiently 
small seed BHs in the centres of the galaxies, we prefer to 
choose their initial mass, Mseed, thereby introducing a free 
parameter. We set Mseed = 10^ M© in agreement with pre - 
vious cosmological simulations (e.g. iBooth ScW3l2009l ). 
However, BHs are still spawned only in cells belonging to the 
maximum level of refinement and that verify equation ([T]) . 
One consequence of this self- controlled formation of the sink 
BHs is that they are not allowed to accrete gas when the 
Jeans criterion is violated. The only way for them to accrete 
gas is to do so by a reasonable physical process such as Bondi 
accretion. With this prescription for initializing the mass of 
the seed black hole, it is conceivable that gas could be nu- 
merically violently Jeans unstable, but this issue is partially 
solved by the consumption of gas in the star forming process 
that temporarily restores gravitational stability. 

To get only one BH per massive galaxy, a halo finder 
is usually run on-the-fly during the simulation to check if 
candidate galaxies alre ady host a BH (|Di Matteo et al ] l2005l : 
iBooth Schavell2009l ). We prefer a simpler, more direct, 
and computationally cheaper approach. To avoid creating 



multiple sinks inside the same galaxy, we ensure that each 
time a cell could potentially produce a sink (i.e. it verifies 
eq. ([1])), it is farther than a minimum radius rmin from all 
other pre-existing sinks. This distance has to be larger than 
the typical size of galactic discs and smaller than the typical 
average inter-galactic distance. Test runs suggest that the 
choice rmin = 100 kpc produces very satisfactory results. 

To avoid formation of sink particles in low density re- 
gions that are Jeans-unstable, we set a minimum threshold 
for the density p > po of gas that can create a new sink, 
where po is the same density threshold that we use for star 
formation. In order to avoid the creation of sink black holes 
before the formation of the very first stars, we check that 
the star density verifies 



/* 



> 0.25, 



(2) 



p* + p 

before a sink is spawned, where p is the gas density. 

When a sink particle is finally created it is split into 
several cloud particles with equal mass. Cloud particles are 
spread over a 4 Ax radius sphere and positioned every 0.5 Ax 
in (x,y,z). The exact number of cloud particles in this con- 
figuration is therefore ricioud = 2109 per sink. This splitting 
process is useful in many ways. First, it keeps a heavy sink 
particle from becoming the dominant gravitational contri- 
bution in one single cell. The latter could catapult particles 
far from their host galaxy by two body encounters. Second, 
it provides a simple canvas over which to compute averaged 
quantities around the sink, for example we use it to deter- 
mine the Bondi accretion rate. 

Sinks are also allowed to merge together if they lie at a 
distance closer than 4 Ax from each other. Mass is conserved 
in this process and momentum vectors of the old sinks are 
simply added to compute the momentum of the new sink. 

Finally we insist on the fact that sink positions and 
velocities are updated in the classical way used to update 
standard particles such as DM particles. No correction on 
their positions and velocities is done to force them to stay 
near their host galaxy (as could be done with the Halo finder 
approach). Thus, weakly bound BHs, such as BHs in satellite 
galaxies of large groups and clusters, may easily be stripped 
from their host galaxy. These BHs behave like star parti- 
cles that tidal forces compel to populate the stellar halo of 
massive galaxies. 

2.3 Accretion rate onto SMBHs 

Since we fail to resolve the accretion disk around our 
SMBHs, whose size is sub-parsec even for the most mas- 
sive ones (^ 10""^ pc according to iMorgan et alll2010l from 
micro- lensing estimates), we use the common prescription 
that these B Hs accrete gas at a Bondi-Hoyle-Lyttleton rate 
(lBondilll952h 



Mb 



47TaG^Mi^p 
(ci+^;2)3/2 



(3) 



where a is a dimensionless boost factor {a ^ 1), Mbh is the 
black hole mass, p is the average gas density, Cs is the aver- 
age sound speed, and v is the gas velocity relative to the BH 
velocity. One of the major difficulties encountered with the 
computation of the relative gas velocity is that in cosmo- 
logical runs, the ISM is poorly resolved and leads to a very 
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thin scale height for galaxies (compared to the resolution). 
Moreover, due to poor sampling of the gravitational force 
in the galactic disc, BHs can slightly oscillates in their host 
galaxy. For this reason a BH close to the centre of a galaxy 
can feel the infalling material coming from the halo at a rel- 
ative velocity much higher than the typical velocity inside 
the bulge. As v is not a reliable quantity, we do not com- 
pute V as the gas velocity relative to the sink velocity but we 
prefer to set it to the average gas velocity dispersion in the 
ISM which is assumed constant and equal to a = lOkm.s"^ 
(|Dib et al.ll2006l : lAgertz et al.ll20Q9l ). 

The average density p and sound speed Cs are computed 
around the BH using the cloud particles for this operation, 
as mentioned in section 12.21 To compute the averages, the 
cell in which each particle sits is assigned a weight ^iven by a 
kernel function w, similar to the one used in lKrumholz et al] 
(I2OO4I ): 



to 



w oc exp 



(4) 



where r is the distance from the cloud particle to the sink 
particle and vk is the radius defined as 




tbh < Ax/4 , 
Ax/4 ^ rsH ^ 2Ax , 
tbh > 2Ax . 



The Bondi-Hoyle radius tbh is given by: 

GMbu 
rsH = 9 — 5 



(5) 



(6) 



where Cs is the exact sound speed in the cell where the sink 
lies. 

The true accretion rate onto the sink is finally limited 
by its Eddington rate 

• 47rGMBHmp 
MEdd — , ( < ) 

CrCTTC 

where ctt is the Thompson cross-section, c is the speed 
of light, mp is the proton mass, and Cr is the ra- 
diative efficiency, assume d to be equal to 0.1 for the 
IShakura SunvaevI (|l973l ) accretion onto a Schwarzschild 
BH. 

The accretion rate is computed at each time step and 
a fraction Mbh At/ricioud of gas mass is depleted from the 
cell where the cloud particle lies and is added to that cloud 
particle, and its sink mass is updated accordingly. At each 
coarse time step cloud particles are re-scattered with equal- 
mass Mbh /'/T'cioud • As the timestep does not depend on the 
accretion speed onto BHs and as low density cells can be 
close to high density cells, a BH might remove more mass 
than is acceptable. To avoid negative densities and numer- 
ical instabilities arising from this, we do not allow a cloud 
particle to deplete more than 25% of the gas content in a 
cell. 

In such large scale cluster simulations, it is impossible to 
resolve the scale and the dumpiness of the ISM. To prevent 
the collapse of the gas from numerical instabilities and to 
take into account the mixing of the different phases in the 
ISM (cold and warm components), we use the polytropic 
EoS described in section 12.11 Applying this EoS means that 
it is impossible to know the "true" density and the "true" 
sound speed in the ISM, thus the accretion rate onto the BHs 
must be modified. Previous work modeling the accretion rate 
onto BHs with such a polytropic EoS set the a parameter 



a constant 100 (ISpringel et aLlliooBI : ISiiacki et aDl2007l : 



Di Matteo et al 



Booth Schave 



I2OO8 ). he prescription from 

(I2OO9I ) who show that a = (p/po)^ is the 
best parametric choice to match observational laws. 

We stress that this polytropic EoS has important conse- 
quences on the accretion rate onto BHs in gas rich systems: 
equation ^ turns into Mbh oc M|hP^^^5 and the temper- 
ature dependence is removed. On the other hand, as soon 
as the cold gas component has been evaporated by star for- 
mation or feedback mechanisms from massive galaxies, the 
accretion rate of the black hole is, by definition, the proper 
Bondi accretion rate. This a boost of the accretion rate is 
an artificial way of modeling the very fast accretion of gas 
within cold and gas-rich galaxies at early epochs, where the 
dumpiness of the ISM is unresolved in large-scale cosmolog- 
ical simulations. 



2.4 AGN feedback modeling 

Followine; lOmma et all (|2004h . we assume that the primary 
source of AGN feedback is a sub-relativistic, momentum- 
imparting bipolar outflow. As discussed in detail by these 
authors, there are many considerations which support this 
hypothesis, the most blatant one being that bipolar out- 
flows are observed around virtually all accreting objects 
in the Universe: stars, black holes and galaxies. Along 
with these authors, we further assume that the advection 
dominated infiow-outfiow sol ution (ADIOS) developed by 
iBlandford BegelmanI (|l999l ) to explain the low luminosi- 
ties of AGNs compared to their estimated Bondi-Hoyle ac- 
cretion rates is correct. The most important feature of the 
ADIOS model, as far as we are concerned, is that the bulk 
of the accretion energy (released as plasma falls into the 
SMBH) drives a (sub-relativistic) wind from the surface of 
the accretion disk. We emphasize that, as pointed out in 
lOmma et af] (|2004l ). this bipolar wind is distinct from ob- 
served relativist ic synchrotron jets which are probably pow- 
ered by the spin of the SMBH itself, although both jets 
can simultaneously be present. However the synchrotron jet 
is very likely irrelevant in terms of AGN feedback since the 
mechanical luminosity of the sub-relativistic outflow is much 
higher than the synchrotron luminosity. Finally, we note 
that in very dense environments it is probable that much of 
the accretion energy is radiated away, driving outflows with 
velocities ~ 0.1c in objects with photon l uminosities also 
on th e order of the Eddington luminosity ([King Poundsl 
I2OO3I ). We argue that as more and more shock-heated rar- 
efied gas fills the central parts of the growing host halo, more 
and more energy will come out as mechanical energy. More 
specifically, we assume that a fraction Cf of the radiated en- 
ergy is imparted to the ambient gas 



^AGN — e^Lr — CferMBHC 



(8) 



in the form of a kinetic jet. 

Such an implementation has the advantage of continu- 
ously releasing energy without radiating it away by cooling 
on the scale of a hydro timestep. Indeed this problem plagues 
supernovae feedback m odelling where energy is generally in- 
jected in thermal form ([Navarro Whit 3l 19931 ) and leads to 
the feedback having no dynamical impact on the surround- 
ing gas. To bypass this issue, some authors release AGN 
thermal energy only after a sufficient amount of gas has 
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Figure 1. Projected DM density (upper panels) and stellar density (bottom panels) in the resimulated 14 Mpc on a side region (left-hand 
panels) and in the cluster, 1.8 Mpc on a side region (right-hand panels). The red circle in the upper right panel shows the rsoo = 940 
kpc radius of the DM halo at z = 0, and green circles in the lower right panel show the virial radii of the stellar structures as detected 
by the halo finder. 



been accreted ont o the BH so as to more severely impact the 
ambient medium ([Siiacki et al. I l2007l :l Booth Sc have "2009^; 
iTevssier et al.ll2010h . In contrast, we model the AGN feed- 
back as a jet-like structure wi th the same momentum profile 
defined in lOmma et al.l (|2004l V This model has already been 
used to follow the self-consistent BH g rowth and its energy 
relea s e in an isolated ga laxy cluster (jCattaneo Tevssierl 
I2OO7I : [Pubois et al]|2009h . Mass, momentum and energy are 
spread over a small cylinder of radius rj and height 2hj 
{hj for one side of the jet) multiplied with a kernel window 
function 



^('■cy.) = ^exp(^-^) ' 

where rcyi is the distance to the axis of the cylinder. The 
mass deposition follows 

Mj (reyi) = ^^vMbh , (10) 

where is the integrated value of ip over the whole cylin- 
der, and ?7 = 100 is an arbitrary value that represents the 
mass loading factor of the jet on unresolved scales. Mass is 
transferred from the central cell (where the sink lies) to all 
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the enclosed cells within the jet. Momentum qj is deposited 
in opposite directions from the centre along the jet axis, 
according to 

||qj|| (rcyi) = Mj (r.„) IIujII = t^M^^,/^c^^ , (n) 

where ||uj|| = (2efer/'/7)^^^c is the velocity of the jet (||uj|| ~ 
lO^km.s"^ for ef = 1), j is the unit spin vector of the BH 
which defines the jet axis, and dr is the distance vector from 
the centre of the black hole, j is computed by adding the 
different contributions from the neighbouring cells (sampled 
with the cloud particles) to the total angular momentum 



J = m^dri X Ui , 



(12) 



where rrii and ui are the mass and velocity of the gas in 
the cell harbouring the cloud particle i, so that j = J/||J||. 
Finally the kinetic energy deposited within a cell is 



Ej (rcyi) 



qj (^cyl) _ 1p (rcyl) 



2Mj (rcyi) 



^AGN . 



(13) 



Integrating this energy deposition over all the cells within 
the jet, we recover Eagn- 

We point out that our jet has no opening angle and 
should therefore propag:ate along a st raigjht line as it is per- 
fectly collimated. lOmma et all (|2004l l have shown that this 
kind of jet stays collimated over quite long distances (100 
kpc) compared to its initial broadness and length (1 kpc). 
It broadens as it reaches equilibrium with the surrounding 
hot a mbient medium. The same beh avior is also pointe d 
out bv lCattaneo TevssieJ (|2007l ) and lPubois et al.1 (|2009l ). 
but with the difference that when strong turbulent mo- 
tions begin to develop in the cluster core due to the for- 
mation of a cooling flow, the jet is more quickly mixed 
with the ICM. The choice of the jet velocity input ||uj|| ^ 
lO^km.s"^ (or equivalently the mass loading factor rj = 
100) is particularly arbi t rary but based on earlier w orks 
from lOmma et al.1 (l2004l \ ICattaneo Tevssierl (l2007l ) and 



iDubois et all (|200' 



The same simulation performed with 
||uj|| ~ 3.10^ km. s~^ produces results in very strong agree- 
ments with the ones presented here (mass of the most mas- 
sive BH and stellar mass of the central galaxy agrees within 
1%) suggesting that even strong variations of ||uj|| keep our 
results unchanged. 

We set rj and hj equal to Ax, and the energy ef- 
ficiency ef = 1 so as to reproduce the Mbh-M* and 
Mbh-ct* observational relations. Larger values of rj and 
hj have been tested at the same resolution Ax c^i 1 kpc 
in a cosmological simulation (as opposed to a resimula- 
tion like the one presented in this paper) and they pro- 
duce BHs which are too massive with respect to their 
host galaxy. Note that our total efficiency CrCf = 0.1 is 
also in good agreement with the average value obtained 
by general relativistic magneto-hydrodynamics numerical 
simulations of the accretion- eject io n mechanism in accre- 
tion discs around spinn ing BHs fe.g.lDe Villiers et"alll2005l . 
iHawlev KrolikI l2006l . or lBenson Babulll20od " and refer- 
ences therein). Lower ef values again cause black holes to 
become more massive, overshooting the Mbh-M* observa- 
tional relation. 



3 SIMULATION SET-UP 

The simulations are run wi th the Adaptiv e Mesh Refine- 
ment (AMR) code RAMSES (|Tevssierl I2OO2I ) . The evolution 
of the gas is followed using a second-order unsplit Godunov 
scheme for the Euler equations. The Riemann solver used to 
compute the fiux at a cell interface is the acoustic solver us- 
ing a first-order MinMod Total variation diminishing scheme 
to reconstruct the interpolated variables from their cell- 
centered values. Collisonless particles (dark matter, stars 
and sink particles) are evolved using a particle-mesh solver 
with Cloud-in-Cell (CIC) interpolation. 

We assume a fiat ACDM cosmology with total matter 
density — 0.3, baryon density = 0.045, dark en- 
ergy density Qa = 0.7, fiuctuation amplitude at 8/i~^.Mpc 
(78 = 0.90 and Hubble constant Hq — 70 km.s~^.Mpc~^ 
that corresponds to the Wilkinson Microwav e Anisotropics 
Probe (WMAP) 1 year best-fitting cosmology (jSpergel et al.l 
j2003j)- The simulations are performed using a resimulation 
(zoom) technique: the coarse region is a 128^ grid with 
Mdm = 2.9 X 10^° M© DM resolution in a 80 h ^Mpc simula- 
tion box. This region contains a smaller 256^ equivalent grid 
in a sphere of radius 20h~^Mpc with Mdm = 3.6 x 10^ M© 
DM resolution, which in turn encloses the final high resolu- 
tion sphere with radius 6h~^Mpc, 512^ equivalent grid and 
Mdm = 4.5 x 10^ M© DM resolution. Fi gure [1] shows the 
distribution of DM and the distribution of stars in the zoom 
region and in the galaxy cluster at z = 0. 

The smallest region is the resimulation zone where cells 
may be refined up to ^max = 16 levels of refinement, reaching 
1.19 h~^. kpc, following a quasi-Lagragian criterion: if more 
than 8 dark matter particles lie in a cell, or if the baryon 
mass exceeds 8 times the initial dark matter mass resolution, 
the cell is refined. This strategy allows AMR codes, such as 
RAMSES, which use CIC interpolation in their gravity solver, 
to avoid prop agating discreetness noise from small scales 
(|Romeo et al.ll2008i V A Jeans length criterion is also added 
to ensure the numerical stabili ty of the scheme on all levels 
^ < ^max (jTruelove et al.lll997l l, and where 5p = p/p > 10^: 
the cells fulfilling these conditions must sample the local 
Jeans length with more than 4 cells. We point out that the 
^max = 16 level of refinement is only reached at aexp = (1 + 
z)~^ = 0.8, and that the actual maximum level of refinement 
for a given redshift is increased as the expansion factor grows 
with time, i.e. ^max — 15 at aexp — 0.4, ^max — 14 at aexp — 
0.2, etc. This allows us to resolve the smallest scales with a 
roughly constant physical size (0.95 < Ax < 1.9h-^kpc), 
rather than a constant comoving size. 

The resimulated region tracks the formation of a galaxy 
cluster with a 1:1 major merger occurring at z = 0.8 . In- 
deed throughout its formation, the cluster chosen for res- 
imulation passes through different dynamical stages: both a 
morphologically perturbed epoch occurring at half the age 
of the Universe, and a relaxed state at late times, which 
permits us to study the different associated states of the 
BH self-regulated growth. Figure [2] shows the dark matter 
halo merger tree history for this cluster. Halos and sub-halos 
are identified and followed using the M ost massive Sub-nod e 
Method (MSM) algorithm described in lTweed et"all (|2009l ). 
which isolates bound substructures. The cluster experiences 
a major halo merger at z ^ 1.7 and two proto-clusters pro- 
genitors merge together earlier at z :^ 3.1. These mergers 
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Figure 2. Dark matter halo merger history of the resimulated 
galaxy cluster. Only the 8 most massive branches and sub- 
branches are shown. Halos are designated with light blue circles 
and sub-halos with dark blue squares. The final merger between 
a sub-halo and its host halo is represented by a solid line connect- 
ing both objects. Note that the cluster main branch 1 experiences 
two major mergers, one with branch 2 at z ~ 3.1 and one with 
branch 3 at z ~ 1.7. However, these mergers are completed at 
z 1.7 and z ~ 0.8 respectively, epochs which coincide with the 
merger of the central galaxies hosted in these halos. Central BHs 
merge later (see Fig [3]). 



end at z ^ 0.8 and z ^ 1.7 respectively, when the cen- 
tral galaxies hosted in the (sub) halos merge. Central BHs 
hosted by these galaxies merge later, as shown in Fig [3] At 
later times, most of the mass growth of the cluster occurs 
through diffuse accretion or minor mergers. In the following, 
we discuss how such events might trigger or halt the AGN 
activity of the central (most massive) BH. 



4 GROWTH OF SMBHS AND THEIR 
ACTIVITY 

The growth of a BH is tightly linked to the accretion his- 
tor y of its host halo ( c.f. c oeval growth scenari o advocated 
bv iMiller enil l2006l and iHopkins et "aD l2007l ) . In princi- 
ple, cold gas which flows directly into the central nucleus 
in a free-fall time will very efficiently grow BHs. How- 
ever this rapid growth might be substantially reduced by 
AGN activity that could expel both energy and material in 
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Lookback time (Gyr) 

Figure 3. Mass evolution (upper panel) and accretion rate rela- 
tive to the Eddington limit (bottom panel) of the most massive 
BH (black solid lines), and two of its most massive progenitors 
involved in the episodes illustrated in figure [H and O i.e. central 
BH of branch 2 (blue dot-dashed lines) and branch 3 (red dashed 
lines) of the merger tree (figure O. Black vertical dotted lines in 
the upper panel mark the mergers of the progenitors with the 
most massive BH. 



the vicinity of the BH. In this case, self-regulation of the 
BH growth possibly drives the relations o bserved between 
BH mass and their host ga l axy properties (iMagorrian et aP 
Il998l : iTremaine et"alll2002l : iHaring Rixll2004l ). 

At high redshift (z > 2), most galaxies seem to har- 
bor a massive cold gas disc compon ent, both in cosmolog- 
ical hydrodynam ics simulations (Oc virk et al.i r2008) and in 
the observations ([Shapiro et ahlboOSl ), which can be tapped 
to fuel rapid BH growth. Accordingly, in the cosmological 
re-simulation of a cluster presented here, the central BH 
reaches a few tenths of its final mass when the Universe 
is less than 2 Gyr old, accreting at a rate above ^1 % 
of its Eddington limit (figure [3|) . As the initial seed mass 
of the black hole is 10^ M0, this means that its mass in- 
creases by a factor 10^ in a tenth of the Hubble time. On 
Fig. [3] the growth of two of its most massive BH progen- 
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itors (indicated by red dashed and blue dot-dashed hnes 
on the figure) is also displayed until they merge with the 
final BH (Mbh = 1.7 10^° M© at z = 0). These merg- 
ers (vertical dotted lines in upper panel of Fig [3]) coin- 
cide with two important halo mergers in the history of the 
cluster (branches labelled 2 and 3 on Fig. [2|). At high red- 
shift {z > 2), these BH progenitors behave like the main 
one: they accrete gas at a fraction of their Eddington rate 
and this fraction stea dily decreases wit h time f rom z = 4 
to z = 2. We know (iBirnboim Deke l 2003; Keres et al.1 
l2005l : iDekel k BirnboimI l2006l : lOcvirk et al.. 120081 ) that for 
the most massive halos (those with masses A^stream ^ 6 X 
10'^(l + 2;)^ Mo at z ^ 2), cold accretion of gas from the IGM 
is efficiently thermalized at a few virial radii by an accretion 
shock, and as a result we expect the accretion rate in the 
centre of the halo to drop. From Fig. [2] one can see that the 
DM halos of branches 2 and 3 have masses ~ 10^^ Mo at 
redshift z = 3 and therefore satisfy the Mstream criterion. 
On that account we claim that this explains, in part, the 
decrease of the BH accretion rate relative to its Eddington 
rate. 

Mergers are a non-negligible growing mode at interme- 
diate and lower redshift, as there is less cold gas to feed the 
BH in the massive cluster. Indeed fig. [3] shows that the BH 
doubles its mass at redshift ~ 1.6 when two BHs of com- 
parable mass (Mbh = 3.1 10^ Mo and Mbh = 8.2 10^ Mo) 
coalesce. The extra amount of mass f:^ 210^ Mo) comes 
from the fast accretion of material brought in by the galaxy 
major merger (M* = 1.3 10^^ Mo and M* = 6.1 10^^ Mo). 
This merger appears in the DM merger tree (figure [2|) when 
branches 1 and 2 join at redshift 1.6. It is interesting to note 
that the first BH that forms is not necessarily the most mas- 
sive one at late times (in our case the BH h osted by branch 2 
halo f orms first), as already pointed out bv lDi Matteo et aP 
(l2008l ). 

Fig. 51 [5] and [6] show three different episodes of the 
formation and evolution of the cluster, respectively a high- 
redshift major merger between two gas-rich galaxies (Fig.[4|), 
the major merger of the two clusters (Fig.[5|), and the relax- 
ation of the cluster at late times (Fig[6|). 

The halo merger between branches 1 and 2 at z = 3.1 
results in a cataclysmic episode for its host galaxies at 2; :^ 
1.6: a large amount of gas is expelled far from the core of 
the halo, reaching the virial radius, and the resulting disc 
of gas from the two merging galaxies is almost completely 
disrupted. This sequence of events is illustrated in figure UJ 
On the left panel, we observe the encounter of the two gas- 
rich galaxies before they merge. Shortly after they merge 
(fig.|3]middle panel), their respective BHs do so as well which 
results in a strong jet that disrupts most of the cold baryon 
content in the galaxy and shock heats the ambient medium 
to high temperature. The jet propagates supersonically at 
Mach 3 (iXjet ^ 3000 km.s"^ and Cs ^ 1000 km.s"^) before 
being stopped by the intergalactic medium at r ^ 1.2 Mpc, 
which corresponds to about 3rvir at this redshift (fig fright 
panel) . 

The disruption of cold material by AGN feedback has 
already been noted by iDi Matteo et al.) (| 20051 ) in idealized 
simulations of a gas rich merger. It is comf orting to confirm 
their results within a cosmological setting. iKhalatvan et al.] 
(|2008l l have also pointed out that mergers could trigger a 
high level of AGN activity during the formation of a small 



galaxy group. Finally, we see two hotspots during the jet 
propagation, that look like radio lobes (fig. |4] right panel). 
Such events (strong jets following a merger) become rarer 
as time goes on since the combined action of star formation 
and early AGN activity strongly diminishes the cold and 
dense gas content in massive halos. 

To check how this striking result is affected by our lim- 
ited spatial resolution, we performed the same simulation 
with one more level of refinement (Ax = 0.6h~^.kpc). The 
same burst appears at the same redshift but its power is 
slightly lower, because more gas has been pre-heated by 
a strong AGN activity in a previous merger taking place 
at higher redshift ^ ^ 4. This is not a very surprising ef- 
fect: with more resolution, the density contrast is more pro- 
nounced especially in poorly resolved high-redshift galaxies, 
which, in turn, leads to a faster accretion rate at early times. 
As in our standard run, the cold gas in the core of the halo 
is strongly disrupted by the AGN activity triggered by the 
wet merger occurring at z = 1.6, so that BH growth and 
AGN luminosity are suppressed for 3 Gyrs (see figure [71 from 
z — 1.6 to z — 0.6). Therefore we can reasonably claim that 
most of the important features describing the BH growth 
are already captured by our standard resolution run. 

Figure [71 also shows the energy that would be released 
by supernova feedback if it were implemented in our simula- 
tion. To estimate this energy release, w e assume that stars 
are distributed according to a lSalpeterl (17955 ) Initial Mass 
Function, for which each massive star (M^ > 8 Mo) deposits 
10^^ erg per 10 Mo into the ISM. We see that the energy 
from this form of feedback is always lower than that from 
the AGN at all times, except during the post-merger phase 
from z = 1.6 to z — 0.6. However, in this post-merger phase 
the high level of supernova feedback is an artifact of the way 
star formation is computed: we trace back the star forma- 
tion history (SFH) of the central galaxy using all the stars 
which belong to it at z = 0, so that the star formation and 
hence the supernova rate we derive from it includes that of 
its accreted satellites. In this redshift range (0.6 < 2; < 1.6), 
the star formation rate is dominated by the galaxy progen- 
itor that has not undergone the cataclysmic quasar phase 
(hosted by the branch 3 halo on Fig [2]), whilst the star for- 
mation activity in the quasar galaxy progenitor (hosted by 
the branch 1 halo) is completely suppressed. In light of this, 
it is a fair approximation to neglect the feedback from su- 
pernovae on the evolution of this galaxy cluster. 

At z ^ 1.7, another halo major merger (1:1) occurs 
(branches 1 and 3 in the merger tree (figure [2])). However, 
the most massive (Mbh = 5.9 10^ Mo) BH only merges with 
the (Mbh = 4.3 10^ Mo) BH of its cluster companion at a 
much later epoch {z = 0.58 right panel of Fig [5]). Note that, 
in this case, the BH merger also takes place quite a long 
time after the central galaxies hosting the BHs merge to- 
gether. As a matter of fact, the 2; = 0.8 major galaxy merger 
drives the cluster gas to temperatures twice the virial tem- 
perature thanks to a violent shock wave (middle panel of 
Fig [5]). However, during this galaxy merger phase, the ac- 
cretion rate onto the most massive BH drops to negligible 
values (^ 10~^MEdd fig- [3 solid black line), and only its 
companion continues to accrete at a moderate rate (a few 
10"^MEdd, fig. [31 red dashed line). The resulting AGN ac- 
tivity, even when boosted by the final BH merger at 2; = 0.58 
does not seriously impact the highly pressurized ICM gas 
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Figure 4. Projected temperature in linear scale along the z direction at z = 1.76, z = 1.56 and z = 1.49 from left to right of the most 
massive cluster progenitors. Branch 1 of the merger tree corresponds to the red circle and branch 2 to the blue circle (see Fig [2]). This 
redshift sequence corresponds to pre-central galaxy merger, central galaxy and BH mergers and post-central galaxy merger respectively. 
The red saturated region corresponds to a temperature T ~ 2.5 — 3 keV. The size of the images is 3.6 Mpc in comoving units. 




Figure 5. Projected temperature in linear scale along the x direction at z = 0.88, z = 0.74 and z = 0.58 from left to right of the 
cluster during its merging phase. Branch 1 of the merger tree corresponds to the red circle and branch 3 to the blue circle (see Fig [2]). 
This redshift sequence corresponds to pre-central galaxy merger, slightly post-central galaxy merger and BH merger respectively. The red 
saturated region corresponds to a temperature T ~ 4 — 6 keV. The size of the images is 3.6 Mpc in comoving units. 




Figure 6. Projected temperature in linear scale along the z direction at z = 0.09, z = 0.04 and z = from left to right of the cluster 
during its relaxed phase. The size of the images is 1.8 Mpc in comoving units. 
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Figure 7. AGN luminosity of the most massive BH at 2; = as a 
function of time (red curve). For comparison, a simple estimate of 
the contribution from supernova activity based on star formation 
rate is given (black curve). Note however that this contribution 
is negligible compared to that of the AGN and is not included in 
the simulation. 
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Figure 8. Volume- weighted cooling time for the run without 
AGN (black) at z = and with AGN (red) at z = (solid), 
z = 0.04 (dashed) and z = 0.09 (dot-dashed). The horizontal 
dotted lines are the Hubble time and the time since the major 
merger of the two central galaxies at z ~ 0.8. 



which confines the jet energy to the very central parts of the 
cluster (right panel of Fig [5]). After the major BH merger, 
the accretion rate onto the central BH becomes extremely 
small (10~^MEdd) due to the complete evaporation of left- 
over cold material by the final outburst of AGN activity. 
Subsequently the BH stays in this almost-dead phase for 2 
Gyr. 

It is striking that the less massive BH progenitor of the 
latter BH merger (at z — 0.58) accretes gas at 0.01 MEdd 
before the merger, whereas accretion onto the most massive 
BH is negligible. The reason for this different behavior can 
be understood by looking at the temperature maps of both 
cluster progenitors (left panel of figure [5]): the most massive 
cluster (branch 1, red circle) is slightly warmer than its com- 
panion (branch 3, blue circle), as it has been pre- heated by 
important quasar activity at earlier redshifts (z^ 1.56). On 
the other hand, the less massive progenitor did not experi- 
ence such strong pre-heating, and as a result has a lower gas 
temperature, and therefore a higher accretion rate during 
the pre-merger phase. 

Finally, the cluster relaxes and the inner halo cold gas 
reservoir gets replenished, fueling a faint accretion onto the 
BH. This translates into low AGN activity. Episodically, 
stronger jets are produced by the AGN in this phase which 
yield small perturbations of the ICM temperature in the 
form of sound waves (figure [6|. These jets are roughly sonic 
(ujet — 1000 km/s and Cs — 1300 km. s~^), and are efficiently 
thermalized by the ICM. As a result they do not propagate 
farther than a few 10 kpc. 

Nevertheless, the part of the jet energy which is carried 
by these sound waves limits the cooling flow in the cluster 
core. However the cooling time in the core (< 100 kpc) is ex- 
tremely short (smaller than a Gyr, see figure [8]), so the cool- 



ing flow eventually develops again and feeds the BH afresh. 
As a result, the BH accretion rate increases from 10~^ MEdd 
to a few 10~^ MEdd at 0, giving rise to late-time AGN 
activity. 

Figure [9] shows different physical properties in a slice of 
gas cut through the AGN jet at z = 0. They show that the 
jet is at low-density and high temperature in good agree- 
ment with high resolution simulations o f jet- format ion (see 
iHeinz et al.ll2006l : ISimionescu et al. I I2OO9I for example) . A re- 
markable feature is the formation of two under-dense cavi- 
ties filled with sonic-jet material. We have computed a simu- 
lated X-ray map of these cavities in 3 different temperature 
bands (see figure [T0|) . These c avities are r eminiscent of the 
ones observed in Perseus A bv lFabian et al. (2006) in which 
a strong cooling core is also visible. In our simulation, an 
extended cooling-core (a few 10 kpc across) is absent, but 
the cooling flow which gives rise to late-time AG N activity 
is clearly present. As in the iFabian et all (|2006l ) observa- 
tions, we interpret the ripples induced by our jet modelling 
as sound waves. This can be seen in figure [11] where ra- 
dial velocities are always sub-sonic. These sou nd waves, pro- 
vided one can dissipate them viscously (Fabi an et al.l 120031 : 
Ruszkows ki et al. l|2004a") can reheat the ICM at distances 
larger than the scalelength of the jet. In our simulation, 
no explicit viscosity is included, but these spherical sound 
waves do not appear at radii larger than rsoopc? suggesting 
that they have been dissipated by numerical viscosity on 
these scales. 
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Figure 9. Density (uppermost panel), temperature (second from 
top), entropy (third from top) and pressure (bottom panel) slices 
through the central AGN at z = 0. Color bars on the side of each 
panel indicate units in H.cm~^, keV, keV.cm"-^, and erg.cm""^ 
from top to bottom. Images are 447 kpc on a side. The cone 
shows the solid angle used to compute the averaged quantities in 
figure 111! 




Figure 10. Composite RGB image of the simulated Xray emis- 
sion 0.3-1 (red), 1-3.5 (green), 3.5-10 keV (blue) bands for the 
AGN run at z = 0. The image size is 447 kpc on a side. 



5 REGULATION OF THE COOLING 
CATASTROPHE 

In the absence of strong feedback processes to offset the 
cooling of gas in the potential wells of massive DM ha- 
los, too many massive galaxies are formed both in CDM 
cosmological hydrodynamical simulations and semi-analytic 
models of galaxy formation and evolution. In our cluster 
zoom simulation, when no AGN feedback is considered, the 
final mass of stars in the central cD galaxy is very high, 
M* 1.710^^ Mo, for a M500 = 2.4 10^^ M© (M200 = 
2.91O^'^M0) dark matter halo with radius rsoo = 940 kpc 
(resp. r2oo = 1370 kpc, see figure [T]l0. In the presence of 
stirring from AGN feedback, the total stellar mass is re- 
duced to M* 5.6 10^^ Mo, i.e. by more than a factor of 
3. To compute t he stellar mass co ntent, we use the same 
MSM algorithm (iTwe ed et alll2009l ) as for the dark matter 
but with different parameters, since stars are more clustered 
than DM particles. This tool efficiently separates one galaxy 
from another, especially the central galaxy from its satellite 
galaxies (see figure [1]). However, the algorithm used in this 
method attributes all the stars present in the dark matter 
halo and not part of satellite galaxies to the central one. As 
a result, a non- negligible part of the stellar mass of the cen- 
tral galaxy resides in the intra-cluster stellar halo (composed 
of all the stars stripped from satellite galaxies of the clus- 
ter), which has a very large extent (up to ^ 400 kpc). This 
caveat must be borne in mind when comparing the stellar 
mass of the central object with observations: our estimate 
only provides an upper limit of the stellar mass content of 
the central galaxy. 

Figure [T2l compares the stellar mass evolution as a func- 
tion of redshift for the most massive galaxy at 2; = (solid 
line) and its two most massive progenitors (dashed and 
dot-dashed lines, central galaxies of DM halos identified as 
branch 1 and 3 in figure [2| in the AGN and no- AGN runs 
(red and black sets of curves respectively). We observe that 

^ All quantities with subscripts 200 or 500 refer to regions 
with overdensities 200 or 500 times larger than critical (pc = 
3H^/(87tG)). 
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Figure 12. Stellar mass evolution of the two galactic progenitors 
(dashed and dot-dashed lines) involved in the major merger of the 
cluster for the run without AGN (black) and with AGN (red). 
The solid lines show the cumulative stellar mass of these two 
progenitors. Crosses and diamonds indicate the 0.5 M* (2; = 0), 
and 0.1 M*(2; = 0) epoch respectively. 
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Figure 11. From top to bottom: volume- weighted average tem- 
perature, pressure, radial velocity and density as a function of 
radius in the solid angle seen in figure [9] at 2; = for the AGN 



Figure 13. Star formation rate as a function of the lookback time 
for the most massive galaxy at 2; = for the no- AGN (black) and 
AGN run (red). 



the reduction of the stellar mass is a continuous process 
which begins at an early stage (around 2; ^ 5) but gets am- 
plified as time goes on to reach a factor 3 at 2; 1. Both 
progenitors seems to follow the same reduction of their stel- 
lar content which suggests that even the branch 3 cluster, 
which does not exhibit any strong quasar activity at high 
redshift, is able to prevent some gas from falling onto the 



central galaxy. The process is more quiescent as can be seen 
in its BH accretion rate (red dashed curve in figure El, but 
even this continuous and moderate AGN activity can effi- 
ciently reduce the star formation. 

We evaluate the bulge to disc mass ratio of the central 
galaxy by using the kinematics of its star particles. First, we 
identify the rotation axis of the galaxy to define the correct 
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Figure 14. Integrated total mass (solid lines), dark matter mass 
(dotted lines), stellar mass (dashed lines) and gas mass (dot- 
dashed lines) for the run without AGN (black) and with AGN 
(red). The dotted vertical line where r = 0.1 x rsoo- 



cylindrical reference frame in which we project the velocity 
components of each star particle. A particle belongs to the 
bulge if its circular velocity is lower than half its total ve- 
locity. With that definition, the bulge-to-total mass ratio of 
the central galaxy is 0.75 for the simulation with AGN, and 
0.80 for the simulation without AGN. Thus it appears that 
although AGN dramatically change the SFH, they have a 
much less significant impact on the morphology of a galaxy. 

Figure [13] shows the star formation rate for the central 
galaxy as a function of time for the two runs. To compute 
the time evolution of the SFR, we simply have identified the 
stars belonging to the central galaxy at 2; = 0, and traced 
them back in time. Thus it is the SFH summed over all 
the stars of all the satellite galaxies that have been accreted 
onto the central galaxy throughout its evolution. The SFR 
continuously decreases with time due to early AGN activ- 
ity (z^ 4), but the dramatic decline in SFR occurs around 
z = 0.6 when the BHs hosted by the central galaxy merger 
remnant of branch 1 and 3 halos finally merge. The vast ma- 
jority of the cold gas is heated up by the AGN during this 
violent merger. In contrast, the cold gas in the no- AGN case 
is simply compressed in the galaxy merger which results in a 
double small star formation peak around the same redshift 
(z ^ 0.6). The latter effect i s a well-known property o f merg- 
ing galaxies without AGN (iMihos HernauistlfToOGl ). How- 
ever, mergers of galaxies containing BHs boost the accretion 
of gas onto the BH fueling strong AGN activity and produce 
a dip in the SFR by reducing the cold-gas content. These 
features are clearly seen at redshift z = 1.6 and z = 0.6 in 
figure 1131 Such a behavior has already been analyzed in de- 
tail in idealized (as opposed to cosmological) simulations o f 
galaxy merger (ISpringel et al.ll2005l : [Di Matteo et al.ll2005l ). 
Our work confirms that it also occurs in more realistic cos- 
mological configurations. 

We have measured the cumulative mass profiles of the 



Figure 15. Baryon (solid lines), dark matter (dotted lines), stel- 
lar (dashed lines) and gas (dot-dashed lines) cumulated mass frac- 
tions for the run without AGN (black) and the run with AGN 
(red). Dark blue is the cold gas component and green the hot 
gas component for the no-AGN run. Light blue is the cold gas 
component and orange the hot gas component for the AGN run. 
The dotted vertical line where r = 0.1 x rsoo- The horizontal line 
indicates the Universal baryon fraction = Q]^/Qrn = 0.15. 



different components (gas, stars, dark matter) as a function 
of radius for the cluster at 2; = ffigure [T¥)) . In the absence of 
feedback, most of the cold baryons are concentrated within 
the galaxy (in the inner 10 kpc): i.e. the cooling catastrophe 
has occurred. There is a strong difference in the cumulative 
profiles between the runs with and without AGN activity, 
especially in the central region of the cluster. They differ by 
a factor 5 in total mass at r = 10 kpc and by a factor 1.5 at 
r — 100 kpc. Without AGN, the gravitational potential is 
steeper in the centre of the cluster, baryons accumulate and 
gravitationally pull DM along with them. This is predicted 
to have severe consequences when sim ulating the gravita- 
tiona l lensing effect of such structures (|Peirani et ahllioosl : 
Meneghetti et al.ll2010l ). 

Baryon fraction as a function of radius provides us with 
a good benchmark to quantify the influence of feedback pro- 
cesses. We have computed the cumulative fractions of stars, 
gas and dark matter for both runs (figure [T5)) . We observe a 
very small difference in the total baryon fraction at r2oo and 
rsoo which means that this quantity is relatively immune to 
the presence of feedback. In contrast, the baryon fraction 
has decreased by a factor 2 at r = O.lrsoo- This means that 
half of the baryons that were concentrated in the inner parts 
of the halo have been redistributed by the AGN feedback at 
larger radii. However, as the baryon fraction seems relatively 
independent of the presence of the AGN at very large radii 
(> T500), we can conclude that only a modest fraction of the 
gas is expelled by the AGN outside of the cluster. 

Stellar and gas fractions yield more clues as to the 
impact of AGN feedback on the cluster history. There is 
a strong decrease in the stellar fraction, even at large 
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Table 1. Comparison of the stellar, gas and baryon fractions at 
r5oo in our simulations with the observational data. 





no AGN 


AGN 


L03« 


G07^ 


G09^ 


fstar 
•^500 

fgas 

•^500 
•^500 


0.090 
0.073 
0.163 


0.036 
0.116 
0.152 


0.019 
0.117 
0.136 


0.023 
0.101 
0.124 


0.032 
0.090 
0.122 



" Observational data from Lin et al.| J2003 

^ Observational data from ..Gonzalez et al.l (|2007|). Their gas frac- 
tion is the best fit to data from IVikhlinin et aP (|2006l ) and 
iGastaldello et al.l (|2007l ). 

^ Observational data from iGiodinil (l2009l V Their stellar frac- 
tion are t he be st fit of their data combined with data from 
iLin et ahl (|2003|). Their g a s fraction is t he be st fi t to data 
from l Vikhlinin et al.l (|2006l V [Arnaud et al.l (|2007l ) and ISun et al.l 
(l2009f) 



radii (r > rsoo), meaning that star formation has been 
efficiently suppressed by the AGN activity. We see that 
the stellar fraction at rsoo has been lowered by more 
than a fa ctor 2, which is co mparable with the results ob- 
tained bv lDuffv et al ] (|2010l ). Tabled show the stellar frac- 
tions /loo^, gas fractions /|qq and baryon fractions (gas 
and stars) /500 measured at rsoo in our simulations and 
compared to X-ra y or near infrared measurements made 
by various groups (iLin et al. I2OO3I: Gastaldello et al. l2007l: 



by various groups (ILiin et al.l IzUUdI: lUastalQello et al.l IzUU /I: 
Vikhhnin etaD I2OO6I: lArnaud et all l2007l : I Gonzalez et al.l 



20071 : iGiodinillioogI : ISun et aLllioOol ) . Observational data val- 



ues are obtained using the best fit these authors provide for 
floQ 5 /foo /500 as a function of M500 • From this compari- 
son with observations, the simulation which does not include 
AGN feedback is clearly ruled out, as the stellar fraction is 
too high by a factor ^3 — 4 and the gas fraction too low by 
25% in the most favorable case. This clearly indicates that 
the ICM has undergone a cooling catastrophe: too much gas 
has been depleted and transformed into stars. On the other 
hand, the simulation with AGN feedback shows reasonable 
agreement with the stellar fraction estima ted from X-ra y 
measurements bv , Gonzalez et al.. ^2007) and lGiodinil (|2009l ), 
but still overestimates t he stellar fraction from near infrared 
data bv lLin et al.1 (|2003l ) by about a factor 2. Gas fraction in 
the simulation with AGN, is also within the range of values 
inferred from X-ray gas emission, albeit on the high side. 
The major flaw of our simulations is their inability to match 
the lower than Universal baryon fraction observed in small 
galaxy clusters (M500 < 10^^ M©). We note that the dis- 
crepancy would be even more blatant h ad we used WMAP 
5 year parameters (|Dunklev et al.ll2009l ) since the Universal 
baryon fraction goes up to 18% for this cosmology. In the 
case of our AGN simulation, the cluster baryon fraction is 
close to the Universal baryon fraction Qb/^m = 0.15, but we 
fail to push gas far enough out of the cluster potential well 
to match lower observational values sitting around 0.13. 

The gas fraction behavior is somewhat counter- 
intuitive: it is larger in the run with AGN feedback, what- 
ever the radius is. This is explained by the fact that AGN 
removes gas from the central regions of the cluster to replen- 
ish its outer parts. AGN feedback thus transforms cold gas 
contained in the central galactic disc into hot and diffuse 
halo gas. Moreover, the gas fraction has a remarkable fea- 
ture in the form of a pronounced dip at intermediate radius 



(15 kpc for the no-AGN and 100 kpc for the AGN case) 
which marks the transition between the cold/dense phase 
(n > 0.1 cm~^), and the hot/diffuse component. Such dips in 
the g as fraction also are co mmonplace in X-ray cluster sur- 
veys (|Vikhlinin et al.ll2006l ). Another interesting result from 
Fig. [15] is that the dark matter to total mass fraction in the 
cluster core (r < lOkpc) is higher in the case of AGN feed- 
back, even though the total amount of dark matter is lower 
in this case. We attribute this to the domination of the mass 
budget by the stellar component which pulls DM along wi th 
it through adiabatic contraction (|Blumenthal et al]|l986l V 

Finally, we compare the thermodynamical properties of 
the gas in the two runs in flgure [TBI We have fltted the gas 
density proflle in the AGN run with a /3 proflle of the form 



: (1 + {r/r,f 



-3/3/2 



(14) 



where ps = 0.5 cm~^, Tc = 10 kpc and /3 = 0.6. This proflle 
matches the density proflle of the relaxed cluster at different 
times {z — 0, z — 0.04 and z — 0.09, which are separated by 
500 Myr) in the intermediate radius range 0.05-1 rsoo- When 
the cluster is relaxed, the same analytic proflle extends to 
the core of the cluster, but as soon as a cooling flow devel- 
ops it fails to describe the numerical gas density accurately. 
Indeed, the core of the cluster shrinks as gas flows in, so Vc 
drops, but the (5 index stays identical as the proflle remains 
unaffected by the cooling flow on large scales. By contrast, 
fltting the gas density proflle with a f5 proflle for the simula- 
tion where no AGN is present turns out to be an impossible 
task, since the core radius becomes smaller than the spatial 
resolution in that case. 

Surprisingly, the simulation without the AGN, which 
has endured a cooling catastrophe for Gyrs, exhibits a very 
hot gas core (temperature in excess of 9 keV) with a very 
steep proflle (see second panel from the top on Fig |16p . Actu- 
ally a massive central cold disc component is also present in 
this run and would appear on Fig [16] if we were measuring 
mass- weighted instead of volume- weighted quantities, but 
the properties of the gas at the center of the cluster would 
then reflect the properties of the ISM instead of the diffuse 
ICM. This very hot thermal part, in the no-AGN run, arises 
from the cluster need for more thermal energy to support the 
gas against the extra gravitational compression generated by 
adiabatic contraction. Due to the very high temperature and 
a lack of diffuse gas around the post-merger galaxy this en- 
ergy is not easily radiated away (the cooling time is greater 
than 2 Gyr as shown on Fig [8]). However we can clearly see 
in flgure [TT] that the gas, close to the galactic disc (r < 50 
kpc), still collapses onto that galaxy due to the lack of pres- 
sure support (flgure [16]), which explains the depletion of the 
gas component at r :^ 100 kpc. 

On the other hand, when the AGN is active the tem- 
perature proflle is stabilized and looks quasi-isothermal in 
the range 0.05-1 rsoo (second panel of Fig [16]). Before 2; = 0, 
the temperature is a factor 2 to 4 higher in the inner 10 
kpc, due to heating from the jet, which remains conflned 
in that region. As the gas radiates away the jet energy, its 
temperature drops and a cooling flow develops because of a 
lack of pressure support in the core (second panel from the 
bottom on Fig [16]). Small variations of temperature with ra- 
dius in the form of wiggles can be observed in Fig [16] (second 
panel from the top, solid red curve). These correspond to the 
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Figure 16. From top to bottom: volume-weighted average num- 
ber density, temperature, pressure and entropy as a function of 
radius for the run without AGN (black solid) at z = 0, and the 
run with AGN at z = (red solid), z = 0.04 (red dashed) and 
z = 0.09 (red dot-dashed). The blue line in the number density 
plot is a /3 profile fit with (3 = 0.6. The dashed black line in the 
entropy plot describes a K oc r^-'^ power-law. 



propagation of sound waves into the intra-cluster medium. 
These waves contribute to reheating the cooling plasma in 
the cluster as a whole by propagating and isotropising the 
energy injected by the AGN jet. They manage to offset the 
extremely short cooling time within r < 0.1 rsoo which is at 
least one order of magnitude shorter than the time elapsed 
since the last major merger (see figure (8]), and thus prevent 
most of the gas from collapsing onto the central galaxy. 

The pressure cavities seen on figure [9] and in the X- 
ray map (figure llOp , are visible in the pressure profile of 
figure [16] (second panel from the bottom): at z = 0, there is 
a small depression in the pressure profile around r ^ 15-30 
kpc, which does not appear at earlier times when the AGN 
is not active enough to form these cavities. This feature is 
also detectable in the radial velocity profile of the gas on 
figure \l7\ (bottom panel) : there is a net radial gas outflow 
at r ^ 15 — 30 kpc whose maximum corresponds to the 
maximum extent of the jet and whose outwardly decreasing 
profile reflects the pre-shocked cocoon region. The volume- 
averaged velocity which we plot on this figure is however 
under-estimated, because the sonic outflowing component 
of the jet is mixed with the quasi-steady flow or inflowing 
regions. The velocity inside the jet is much faster, around 
1000 km/s. 

Finally, entropy profiles (bottom panel of Fig [16]) shows 
a plateau in the cluster core (at r < 20 kpc for the AGN 
run, and r < 300 kpc for the no-AGN run) with a strong 
departure from the scaling power law K oc r^'^. This indi- 
cates the level of turbulent mix ing in the gas as explained in 
detail by iMitchell etHI (|2009l ) and is nicely illustrated by 
comparing entropy profiles on figure [16] (bottom panel) with 
radial velocity dispersion profiles on figure [17] (top panel) . 
Such a comparison clearly shows that the stronger the tur- 
bulence level (or equivalently the radial velocity dispersion) , 
the higher the entropy. 



6 DISCUSSION AND CONCLUSION 

We have numerically studied how a self-regulating model of 
supermassive black hole growth and AGN feedback impacts 
the formation history of a large cluster of galaxies. Using a 
resimulation technique to explicitly account for the cosmo- 
logical context (ACDM Universe) which drives the cluster 
growth and an AMR technique to solve the equations of hy- 
drodynamics without running into entropy issues, we find 
that: 

- BHs accrete gas efficiently at high redshift (z > 2), sig- 
nificantly pre-heating proto-cluster halos in the process. 

- some, but not all, wet (gas-rich) mergers fuel strong 
episodic jet activity which transport gas from the cluster 
core to its intermediate/outer regions. 

- reduced infall of cold gas during the more secular phase 
evolution of clusters produces smaller outbursts from the 
central AGN which contribute to heat the whole cluster via 
sound waves but are inefficient at redistributing the gas out- 
wards. 

- late-time AGN activity forms two large cavities corre- 
lated with the emergence of a small cooling flow. 

Whilst our model for accretion onto the black 
hole is commonly used in simulations in the literature 
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Figure 17. Same as figure[T6]for the volume- weighted average ra- 
dial velocity dispersion (upper panel) and radial velocity (bottom 
panel) of the gas. 



(ISpringel et~all l20Q5l: IPi Matteo et alJ l2QQ5l: ISiiacki et aD 



20071: iDi Matteo et alJ l2008l : iBooth Schavd l2009l : 
Tevssier et al.ll201Ql l, this is the first time, to the best of 



i 



our knowledge, that a momentum driven jet is implemented 
as AGN feedback in cosmological simulations and followed 
in a self- consistent way. We argue that this is one step 
in the correct direction since powerful jets are observed 
in the Universe on scales well res olved by any of these 
simulations (e.g. iBridle et al.l Il994h . Other authors have 
adopted a more phenomenological approach where energy 
is either accu mulated by the BH before bein^ released as a 
thermal pulse (|Booth k Schavd l2009l : [Tevssier et alll2010h 
or simply used as a continuo u s heat source ([Springel et al.l 
I2OO5I : iDi Matteo et al.1 l2Q05l. [ 2008h . or used both heat- 
ing modes jsiiacki et al.l I2007T ). Injection of non-thermal 
relativistic protons in rising buoyant AGN bubbles has 
also been explored as an alternative feedback mechanism 
(|Siiacki et al.l l200"8l V It is worth noting the recrudescence 
of efforts done to improve models of AGN feedback in 
idealised (non-cosmological) simulations of galaxy evolu- 
tion, where kinetic energy is deposited e ither isotropically 
(jPebuhr et al.1 |201Q| : IPower et all I2OIOI) or - as in the 
prese nt paper - as a collimated jet (|Navakshin Powerl 
I2OIOI V These models probably capture more of the relevant 
jet physics that what we achieve here, and we believe that 
as cosmological simulation spatial resolution increases, 
more physical insight into the impact of AGN feedback on 
the ICM will be gained by coupling them to such models. 
While we believe that most of the results we get are similar 
to those obtained with the continuous heating model and 



that our feedback seems less efficient at stopping the cooling 
catastrophe than the accumulated heating prescription, the 
devil certainly is in the details and we defer a more detailed 
comparison to a future paper (Dubois et al in prep). It is 
however interesting to note that all these models, including 
ours, are calibrated to provide an acceptable Mbh vs 
(or Mbh vs ct^) relation so that it is very unlikely that this 
relation can be used to constrain feedback mecha nisms. 

From both semi-analy tic prescriptions (fBower et al.l 
I2OO6I : Cattaneo et al.l 2006) and recent numeri cal simula- 
tions jRuszkowski Ohll2010l : IPuffv et al.ll2010l ). it has be- 
come clear that AGN feedback must play a major role in 
reducing the central galaxy stellar mass in groups and clus- 
ters by a large amount. Our implementation of feedback 
succeeds at least partially in reaching this goal: the stellar 
mass of the central object is reduced by more than a factor 
3 in the run where AGN feedback stirs the gas whose prop- 
erties (density profile, temperature, radial velocity) also are 
in better agreement with observations. Another success of 
our AGN modeling is its capacity to reproduce double cav- 
ities separated by a cold component as seen in the X-ray 
emission of observed clusters (e.g. in Perseus, Fa bian et al.l 
I2OO6IV Our experirnents w ith a simple isotropic thermal in- 
put ([Tevssier et al.|[2010l ) suggest that these cannot be re- 
produced. 

One drawback of our simulation is that no supernova 
feedback is included. However we do not expect this feed- 
back to be energetically relevant as a simple estimate based 
on the star formation history of our galaxies shows that it is 
always an order of magnitude lower than the AGN energy 
input (figure [7|), On the other hand, supernovae feedback 
releases metals into the ICM. As the cluster gas is heated 
to very high temperatures (^ 3 keV), cooling in the ICM 
is primarily due to free-free collisions, thus we do not ex- 
pect metals to strongly alter the cooling rate of the plasma. 
Indeed, at T = 3.5 keV, the relative difference in the net 
cooling rate between a zero and a one third of solar metal- 
licity plasma is 0.2 according to iTozzi Norman[ ([200lh . 
This means that, assuming a Z= 0.3 Z© metallicity in the 
ICM, our simulation underestimates the gas cooling rate by 
20 %. 

As always when analysing numerical simulations, one 
has to worry about spatial and mass resolution re- 
lated issues. Indeed, accretion onto a SMBH happens on 
(sub)parsec scales, well below the kpc size of the smallest 
cell in our cluster re-simulation. Even th ough we use the 
empirically well motivated subgrid model o f iBooth Schavd 
(200^) to account for this lack of resolution, one would like 
to test its validity by performing direct simulations. Whilst 
this is beyond the reach of the current generation of super- 
computers, we have tried to assess both the robustness of 
our subgrid and our jet implementations by performing a 
sub-kpc run and conclude that our results are by-and-large 
unchanged at this increased spatial resolution. This is also 
true for a modest increase in mass resolution. We are there- 
fore quite confident that the conclusions we draw in this 
paper are robust vis-a-vis resolution issues and defer a more 
thorough resolution study to future work. 

Finally, other physical mechanisms, which we do not 
model here, could potentially play a significant role in pre- 
venting the development of massive cooling flows. These 
mechanisms involve tapping into the heat reservoir provided 
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by the outer regions of galaxy clusters to raise the gas tem- 
perature in their core. In particular, recent efforts have been 
made to investigate the importance of ani sotropic thermal 
cond u ction in idealized ^a l axv clusters (^Parrish QuataertI 
I2OO8I : IParrish et alj l2009l : iBogdanovic et alJ I2OQ9I ). These 
studies have shown that the HBI can reorient the mag- 
netic field lines in the cluster core in an azimuthal con- 
figuration that stops the inwa rd heat flux. T he most 
recent of these simulations (Parrish et alJ I2OIOI I demon- 
strated that if some small turbulence is brought to break 
this magnetic field topology, heating is able to proceed. 
One has to wonder, in the context of anisotropic ther- 
mal conduction, whether the turbulence induced by small- 
scale and large-scale motions is able to reorder mag- 
netic fields in cosmo l ogical simulations of galaxy clus- 
ters ("Polag et al."2QQ5VAsai et al."2QQ ^jlDubois Tevssierl 
[2008a; Pfrommer & Jonathan Dursi 201^), and, thus, break 
the shelf-shielding of the heat-flux. In particular, the pres- 
ence of AGN stirring can he lp to disturb this magnetic equi- 
librium (iDubois et al.ll2009l ). so we believe that this problem 
needs to be addressed with MHD cosmological simulations 
including both AGN stirring and anisotropic thermal con- 
duction. 
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